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Abstract Plasmaspheric refilling on an L = 4 flux tube was studied by using a time-dependent, 
hydrodynamic plasmaspheric flow model in which the ion streams from the two hemispheres are 
treated as distinct fluids. In the model the continuity, momentum, and energy equations of a two- 
ion (0 + and H + ), quasi-neutral, currendess plasma are solved along a closed geomagnetic field 
line; diffusive equilibrium is not assumed. Collisions between all stream pairs and with neutral 
species are included. The model includes a corotating, tilted dipole magnetic field and neutral 
winds. Ionospheric sources and sinks are accounted for in a self-consistent manner. Electrons are 
assumed to be heated by photoelectrons. The model flux tube extends from 200-km altitude in 
one hemisphere to 200-km altitude in the other hemisphere. Initially, the upwelling streams 
pass through each other practically unimpeded. When the streams approach the boundary in the 
conjugate ionosphere, a shock develops there, which moves upward and dissipates slowly; at 
about the same time a reverse shock develops in the hemisphere of origin, which moves upward. 
After about 1 hour, large shocks develop in each stream near the equator; these shocks move 
toward the equator and downward after crossing the equator. However, these shocks are probably 
artificial, because counterstreaming flows occur in each H + fluid, which the model can only 
handle by creating shocks. 


1. Introduction 

A long-standing question in plasmaspheric refilling has 
been whether the ion streams that come from each hemisphere 
would form an electrostatic shock pair when they meet at the 
equator or whether they would pass through each other [Banks 
et al , 1971; Schulz and Koons, 1972]. Even if the streams 
pass through each other, it is possible that they could excite 
waves that would change the pitch angle of the ions, leading to 
ion trapping [Schulz and Koons, 1972]. In order to model this, 
however, it is necessary to have a model in which the ion 
streams from conjugate hemispheres are distinguished; this is 
called a two-stream model. 

Previous work involving two-stream models has been de- 
scribed by Singh [1988], Rasmussen and Schunk [1988], 
Singh and Torr [1990], Singh [1991], and Singh and Chan 
[1992]. Singh [1988] applied a two-stream version of the 
model described by Singh et al . [1986] to refilling on an L = 6 
field line. In this model, only the continuity and momentum 
equations were solved; the ion and electron temperatures were 
assumed to be constant. In addition, a Boltzmann relation was 
assumed to hold between electron density and electric poten- 
tial: 
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where n e is the electron density, n Q is the electron density at a 
reference altitude, 0 is the electric potential, which is assumed 
to be zero at the reference altitude, and T e is the electron tem- 
perature. This relation is used to calculate the electric field. In 
this study it was assumed that T e - 10 7, in order to reduce com- 
puter time. Singh [1988] found that the plasma streams pass 
through each other without being seriously affected by the 
other stream. Rasmussen and Schunk [1988] also developed a 
two- stream version of the Singh et al. [1986] model and ap- 
plied it to refilling on an L = 4 field line, but with T t = T,-. 
Results were given both with and without H + -H + collisions. 
There was not much difference between these two cases, but the 
streams were reflected at higher altitudes when collisions were 
included. Rasmussen and Schunk [1988] find that the streams 
initially interpenetrate. The density fronts steepen into 
shocks as they travel toward the conjugate ionosphere; these 
shocks pass through the boundaries. Thereafter, reverse 
shocks develop in the hemisphere of origin for each stream, 
travel toward the conjugate hemisphere, and reverse direction. 
The subsequent motions of the shocks give results similar to 
those of one-stream models, with refilling between the 
shocks. They conclude that refilling eventually occurs from 
the equator downward. 

Singh and Torr [1990] described a two-stream model in 
which the temperatures are allowed to be anisotropic. In this 
model the ion continuity, momentum, and parallel and perpen- 
dicular energy equations were solved. However, production and 
loss processes, heat flow, and collisions were neglected. In ad- 
dition, a Boltzmann relation was assumed to hold between 
electron density and electric potential and was used to calculate 
the electric Held. The electron temperature was assumed to be 
constant. Singh and Torr [1990] used this model to study 
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plasmaspheric refilling in the presence of equatorial transverse 
ion heating, on an L = 4 field line. This transverse heating rep- 
resented the effect of wave-particle interactions. The trans- 
verse heating results in a temperature anisotropy with T x > T\\ f 
which can stop the interhemispheric plasma flow for even 
moderate heating levels. This launches a shock wave, which 
propagates from the equator to the ionosphere from which the 
flow originated. Thereafter a quasi-steady state prevails. Tj_ is 
greater than T\\ within 15 degrees of the magnetic equator, but 
at higher latitudes the anisotropy is reversed. The density pro- 
file has a broad maximum near the equator and minima at mid- 
latitudes. There is also a local minimum right at the equator. 
Singh and Torr [1990] claim that the equatorial density maxi- 
mum is analogous to the observations of trapped ions from DE 
1 [Olsen, 1981; Olsen et al , 1987]. However, it is inconsis- 
tent to include large heating levels without allowing for the 
possibility of heat flow; Singh and Chan [1992] argue that the 
classical form of heat flow is not valid, since the mean free 
path is much greater than the scale length for temperature vari- 
ations. 

Singh [1991] compared results for plasmaspheric refilling 
found by using a two-stream model with anisotropic tempera- 
tures to those found when the temperatures are isotropic. In 
this work, Coulomb collisions were included, and there was no 
equatorial heating. The results are similar in an initial period 
(<5 hours), characterized by supersonic flows from the iono- 
spheres and the formation of shocks after the plasma streams 
pass the equator; these shocks move back toward the equator 
and then down to the ionospheres. Thereafter a temperature 
anisotropy develops with T\\ > causing a downward force, 
which prevents plasma behind a shock from following it as it 
crosses the equator. After the shock crosses the equator, the 
enhancement in r M is much greater than that in T±, with the re- 
sult that the downward pressure gradient force exceeds the up- 
ward mirror force, accelerating the shock downward. This 
shock moves down much faster than in the isotropic case. Due 
to the temperature anisotropy, the flux tube remains depleted 
until Coulomb collisions allow the anisotropy to relax. This 
begins at 7.68 hours, when the anisotropy relaxes near the 
ionospheric boundaries, allowing upward subsonic flow from 
the ionosphere. Singh [1991] found that there was minimal re- 
filling during the initial period of supersonic flow; the major- 
ity of the refilling occurs in the following subsonic stage. The 
subsonic flows were found to last for several days. 

Singh and Chan [1992] studied plasmaspheric refilling on 
an L = 4 field line in the presence of equatorial transverse heat- 
ing with a model that included Coulomb collisions. The early 
time results were not significantly different from those shown 
by Singh and Torr [1990]. At later times. Coulomb collisions 
cause the temperature anisotropy to relax, and upward subsonic 
flows develop, similar to what was shown by Singh [1991]. 
Singh and Chan [1992] also give results from a kinetic simula- 
tion that shows that the electrostatic shocks form at the mirror 
points of the equatorially heated ions. 

In this paper, results will be presented for plasmaspheric re- 
filling on an L = 4 field line; these results were found by using 
a two-stream version of the plasmaspheric flow model, which 
we developed [Guiter et al. , 1991]. This is a time-dependent 
hydrodynamic model in which the H + and 0 + continuity, mo- 
mentum, and energy equations and the electron energy equation 
are solved along a closed dipole field line. The temperatures are 
assumed to be isotropic. The model includes the effects of ion- 
ization, charge exchange, recombination, collisions, and heat 


conduction and allows for external heat sources. Collisions 
with neutral species are included in addition to Coulomb colli- 
sions. Production and loss processes and heat flows are in- 
cluded, in contrast to the model used by Singh and coworkers 
[Singh and Torr , 1990; Singh , 1991; Singh and Chan , 1992]. 
Photoelectron heating of thermal electrons and a neutral wind 
profile are assumed. Although neutral winds do not directly 
affect plasmaspheric density, they do affect the ionospheric 
density profile and thus indirectly the plasmaspheric density. 
The ends of the model flux tube are at 200-km altitude, which is 
low enough so that realistic inflow boundary conditions can be 
used. For this work the initial density profile was similar to 
that used by Rasmussen and Schunk [1988]. Improvements 
over previous studies include realistic ionospheric boundary 
conditions and the influence of 0 + streams (primarily as a 
source for H + ) on plasmaspheric refilling. Results are pre- 
sented from only the first hour of the calculation, because arti- 
ficial shocks develop in each H + fluid by that time, due to the 
presence of counterstreaming flows within each fluid. 
However, interesting features develop in the ionospheric re- 
gion, which has not been included in previous two-stream re- 
filling studies. 

2. Model 

For this work the streams of both H + and 0 + from the conju- 
gate ionospheres were treated as distinct fluids, with all flows 
assumed to be field aligned. Ions that originate from the conju- 
gate hemisphere are called outflowing (out of the plasmasphere 
and into the neighboring ionosphere), whereas ions that orig- 
inate in the given hemisphere are called inflowing (into the 
plasmasphere from the neighboring ionosphere). The time- 
dependent continuity, momentum, and energy equations for 
both streams of 0 + and H + ions were solved along a closed 
dipolar field line. However, the electron gas was treated as a 
one-stream fluid. Electron number density and velocity were 
found by assuming the plasma to be quasi-neutral and current- 
less. The electron energy equation was solved in order to ob- 
tain the electron temperature. These equations were presented 
by Guiter and Gombosi [1990]. 

In this model the collision terms were assumed to have the 
following form: 
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In these equations the summations run over all species, includ- 
ing the neutral species; v st is the momentum transfer collision 
frequency of species s with species /, v so is the production fre- 
quency of species s , and L s is the loss frequency of species s\ 
n 0 , u 0 , and T 0 refer to the density, field-aligned velocity, and 
temperature of the parent neutral species. The velocity- 
dependent correction factors *F st and *s, [cf. Tanenbaum , 1967; 
Burgers , 1969; Schunk , 1977] were calculated for the interac- 
tion of a given H + stream with other ion streams and with 
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neutral species. Collisions between all stream pairs and with 
neutral species were included. Collision frequencies for H + -H + 
and 0 + -0 + interactions were taken from Schunk [1983]; other 
collision frequencies and thermal conductivities were taken 
from Raitt et al [1975]. 

Production and loss processes are similar to those in the 
one-stream version of the model. Ions are assumed to be pro- 
duced only in the hemisphere of their origin, but can be lost at 
any point on the field line. 0 + ions are produced by photo- 
ionization and charge exchange of O with H + . H + ions are cre- 
ated by charge exchange of H with 0 + . 0+ is removed by re- 
combination reactions with N 2 and O 2 and by charge exchange 
with H; H + is lost in charge exchange with O. There is an addi- 
tional production and loss mechanism for H + ions which is 
unique to the two-stream model: charge exchange of H + with H. 
This will be a sink for the ions flowing out of the conjugate 
hemisphere and a source for the ions that originate in the 
given hemisphere. The following expression was used for the 
charge exchange frequency [cf. Banks and Kockarts , 1973, part 
A, p. 223]: 

v C£ (ff + ,ff) = 1.0xlO- I0 n(//)Vr(ff + ) + r„ (5) 

The loss rate for the outflowing ions is given by the product of 
this frequency and the density of the outflowing ions; the pro- 
duction rate for the inflowing ions is equal to the loss rate for 
the outflowing ions. 

The neutral atmosphere model used was MSIS-86 [ Hedin , 
1987] for June solstice during solar minimum conditions. A 
neutral wind profile was derived from the vector spherical har- 
monic (VSH) model [Killeen et al , 1987] and depends on the 
hemisphere. The solar zenith angle and 0 + photoproduction 
rate were found for each computational cell; these were calcu- 
lated by using software taken from a conductivity model devel- 
oped by Rasmussen et al [1988]. The neutral atmosphere 
parameters, solar zenith angles, and photoproduction rates are 
updated every 15 min in order to simulate the rotation of the 
flux tube. Heating of thermal electrons by photoelectrons is 
included in the model; however, the work presented in this 
paper was done for nighttime conditions, so there was no elec- 
tron heating except at high altitudes, where a minimum heat- 
ing rate equal to 2.2 x 10' 13 erg cm' 3 s' 1 was imposed in order 
to match plasmaspheric temperature observations due to 
Comfort et al [1985]. 

The model flux tube extends from 200-km altitude in one 
ionosphere to 200-km altitude in the conjugate ionosphere. 
Boundary conditions for the inflowing ions are set at 200-km 
altitude in the ionosphere of origin. Boundary conditions for 
the outflowing H + ions are set at about 400-km altitude, 
whereas boundary conditions for the outflowing O* ions are set 
at about 1500 km. The altitudes where these conditions are set 
were raised from 200 km to avoid numerical problems; how- 
ever, it should be noted that this can result in reflective bound- 
aries, due to the lower neutral densities at these higher alti- 
tudes. The boundary conditions for the inflowing ions are that 
the ions be in chemical and thermal equilibrium with the neu- 
trals and that the field-aligned velocity of the ions is zero. The 
inflowing 0 + boundary density is found by using photo- 
chemical equilibrium during the day but is set to 1.5 x 10 3 cm" 3 
at night. The inflowing H + boundary density is found by as- 
suming charge exchange equilibrium with 0 + . The boundary 
densities for the outflowing ions are found by extrapolating 


from interior points. For the outflowing H + ions the boundary 
temperatures are found by extrapolating from interior points, 
but for the outflowing 0 + ions the boundary temperature is set 
to the temperature in the cell adjacent to the boundary. The 
boundary velocities are set to the velocities in the cell adjacent 
to the boundary if these velocities are downward; otherwise 
they are set to zero. 

The equations are solved by using a combined Godunov 
scheme/Crank-Nicholson method with dimensional splitting, 
as described by Guiter and Gombosi [1990]. The cell size was 
-12.6 km. The model was run in two time steps: a small time 
step at low altitudes and a larger time step at high altitudes, the 
division between low and high altitudes being about 1500 km. 
The low-altitude time step was 6.25 x 10’ 3 s; the ratio of the 
high- and low-altitude time steps varied between 1 and 5. At 
times, however, it was necessary to use a time step as low as 
3.125 x 10' 3 s for all points on the field line. In addition, to 
prevent numerical problems caused by 0 + shocks at high alti- 
tudes, an artificially high (0.1) momentum transfer collision 
frequency for 0 + ions with neutrals was used when the ratio of 
the number density of 0 + to that of electrons was less than 
0.01; these 0 + shocks were also seen in steady state results 
found with the one-stream version of the model, and were not 
due to refilling. 

3. Results 

The initial condition was chosen to simulate a density de- 
pletion. To do this, the initial H + density profile above 2000- 
km altitude was determined by the following formula: 

where A is magnetic latitude and A* is the magnetic latitude that 
corresponds to 2000-km altitude; a minimum density was also 
imposed, equal to 10 -3 times the density of the other H + stream 
in the first cell in from the boundary. This is similar to that 
used by Rasmussen and Schunk [1988]. The initial density pro- 
files for all species are shown in Figure 1 as a function of mag- 
netic latitude. Diumally reproducible results for June solstice 
conditions at midnight on an L = 4 field line from a one-stream 
version of the model were used to determine the 0 + density 
everywhere and the H + density up to 2000-km altitude for each 
stream in the hemisphere of origin. The 0 + density profiles are 
not necessarily realistic at high altitudes because of the artifi- 
cial resistance used in the one-stream model. The electron den- 
sities were determined by using quasi-neutrality. The initial 
temperature profiles were also taken from the one-stream re- 
sults. The initial velocities were assumed to be zero. 

Because of the density depletion, the ion gases expand up- 
ward into the low-density regions. This can be seen in Figure 
2, which shows velocity and density profiles for the southern 
stream of H + ions between 0 and 15 min. It should be noted 
that positive velocities are upward in the northern hemisphere 
but downward in the southern hemisphere, and that positive 
latitudes correspond to the northern hemisphere. The veloci- 
ties are highly supersonic with respect to the ion-acoustic 
speed C Jt where the ion-acoustic speed is defined by 
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Figure 1. Initial density profiles as a function of magnetic 
latitude used in the two-stream calculation on an L = 4 field 
line. The solid line denotes electrons, the dotted line shows 
the H+ southern stream, the dashed line shows the H + northern 
stream, the dot-dashed line corresponds to the 0 + southern 
stream, and the dash-dot-dot-dot line corresponds to the 0 + 
northern stream. 

At 5 min the ion-acoustic Mach number is approximately 3.2. 
As shown by Singh [1988], if the ion-acoustic Mach number is 
greater than 2.3, counterstreaming is expected to occur. The 
stream has crossed the equator by 5 min and by 15 min the den- 
sity front has began to steepen. In the northern hemisphere 
the velocity changes largely as a result of advectioh of the 
velocity fronts. It can also be seen that there is upward flow in 
both the southern and the northern hemispheres. The upward 
flow in the northern hemisphere is a consequence of the initial 
conditions, which specify a low initial density of southern 


stream ions in the northern hemisphere; these ions flow up- 
ward in response to the upward electric field. The density front 
passes into the northern hemisphere practically unimpeded. 

Similar results were found for the stream from the northern 
hemisphere, with the density front in that case passing from 
the northern into the southern hemisphere. The electron den- 
sity fills in rapidly; by 15 min the equatorial electron density 
is approximately 30 cm’ 3 . 

It can also be seen that in this initial period, downward flow 
develops at low altitudes (-400 km). This is a consequence of 
the initial density distribution below 2000-km altitude, which 
corresponds to that for quiet-time nighttime conditions; for 
such conditions, downward flow develops. * 

By 17.5 min the density front has passed through the out- 
flow boundary. Thereafter the density builds up near the bound- 
ary owing to the large downward fluxes. This can be seen in * 

Figure 3, which shows velocity and density profiles at 20 and 
25 min for the southern stream of H + ions. There is also an up- 
ward velocity spike at approximately 860 km in the southern 
hemisphere which is accompanied by a local density deple- 
tion; this develops into a reverse shock, which moves upward 
(a reverse shock is one where the velocity changes from a high 
to low magnitude while the density changes from a low to high 
value). The maximum velocity in this spike is close to the 
sound velocity. 

The development of the velocity spike is shown in Figure 4, 
which shows velocity and density profiles in the southern 
hemisphere between 12.5 and 17.5 min for the southern stream 
of H+ ions; it should be recalled that in this region negative 
velocities are upward. The velocity spike begins to develop as 
early as 12.5 min. The density profile at 17.5 min shows the 
beginning of the density depletion. From this figure it can 
also be seen that in this time period the upward velocity 


H* southern stream 



Figure 2. (a) Velocity and (b) density profiles as a function 
of magnetic latitude for the H+ southern stream. The solid line 
shows the initial profiles, the dotted line gives results after 5 
min, the dashed line gives results after 10 min, and the dot- 
dashed line shows results at 15 min. 
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Figure 3. (a) Velocity and (b) density profiles as a function 
of magnetic latitude for the H + southern stream. The solid line 
shows results at 20 min, and the dotted line gives results at 25 
min. 
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Figure 4. (a) Velocity and (b) density profiles as a function 
of magnetic latitude for the H + southern stream, between *60 
and -40 degrees. The solid line shows results at 12.5 min, the 
dotted line gives results at 15 min, and the dashed line gives 
results at 17.5 min. 

throughout this region increases, driven by the upward electric 
field and by the density gradient force at high enough altitudes. 
Figure 5 shows altitude profiles of electron and ion (southern 
streams) density at 12.5 min in the southern hemisphere. The 
velocity spike develops in the transition region where 0+ is 
the major ion and the H + density increases with altitude, which 
means that the density gradient force is downward and the mo- 
mentum transfer collision frequency of H + with 0 + ions is rela- 
tively small. This is also the region where the electron pres- 
sure gradient changes drastically with altitude, resulting in the 
rapid decrease of the electric field with increasing altitude. In 


time = 12.5 min 



Figure 5. Altitude profiles of electron and ion (southern 
streams) density at 12.5 min in the southern hemisphere. The 
solid line refers to electron density, the dotted line refers to H+ 
density, and the dashed line refers to 0 + density. 


this region the most important forces (for the southern stream 
of H* ions) are the electric, density gradient, and gravitational 
forces; the sum of these forces is upward and significantly 
larger at about 860-km altitude than at higher altitudes. This 
can be seen in Figure 6, which shows altitude profiles of the 
electric, density gradient, and total forces (for the southern 
stream of H + ions) in the southern hemisphere between 15 and 
20 min; it should be noted that positive forces are downward in 
this region and that the force scale for the total force profile is 
different from that for the other two profiles. As a result, the 
velocity increases more at about 860-km altitude. The velocity 
rise causes more flux to leave this region than the ionosphere 
supplies, resulting in a local density depletion. This in turn 
causes the electron density to decrease, which causes the elec- 
tric field to increase as the electron pressure gradient does not 



density gradient force (x 10~ 20 dyne) 



Figure 6, Altitude profiles of the (a) electric, (b) density 
gradient, and (c) total force in the southern hemisphere. The 
solid line shows results at 15 min, the dotted line gives results 
at 17.5 min, and the dashed line gives results at 20 min. 
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change significantly. From Figure 6 it can be seen that the in- 
crease of the electric field with time is partially offset by the 
increase in the downward density gradient force, so that the to- 
tal force increases in only a limited altitude range. The result is 
that the velocity rises more in that altitude range, leading to 
the velocity spike. 

Figure 7 shows altitude profiles of electron and ion 
(southern streams) density at 20 min in the southern hemi- 
sphere. At this time the density profile for the southern stream 
of H + ions has a broad minimum centered at about 740-km 
altitude. This profile is similar to that found by Raitt et al 
[1975] when there is strong upward flow out of the ionosphere, 
as there is in this case. This density structure steepens into a 
shock, which then moves upward. 

By 30 min a shock has developed near the northern bound- 
ary and then moves upward. This can be seen in Figure 8, 
which shows velocity and density profiles at 30 and 40 min for 
the southern stream of H + ions. This develops as the velocity 
is slowed down (starting at about 20 min) in the northern 
hemisphere in the 400- to 500-km altitude range, first by fric- 
tion with neutrals and then by the density gradient force (the 
density increases toward the boundary). The velocity is slowed 
down from a supersonic value, which means that a shock must 
form. By 40 min the shock strength has diminished somewhat, 
although it continues to move upward. Also, by 40 min the re- 
verse shock in the southern hemisphere has developed into a 
forward shock, which propagates with a shock in the northern 
stream of H + ions. Figure 9 shows density results in the south- 
ern hemisphere for electrons and both H + fluids, at 40 and 45 
min; it can be seen that the shocks in the two H + fluids are 
moving together. The motions of these two shocks are coupled 
by the polarization electric field, which at this time is gen- 
erated by the northern stream of H + ions. 

Figure 10 shows velocity and density results between 50 and 
60 min for the southern stream of H + ions. The shock that 
formed near the northern boundary continues to move upward 
and is slowly dissipating; by 60 min it has almost reached 20 
degrees latitude. Also, at 50 min there is upward flow in the 
northern hemisphere of H + ions which come from the southern 
hemisphere. This flow is driven by both the density gradient 
force and the electric field (which is directed toward the equa- 
tor). The shock that originated in the southern hemisphere 


time = 20.0 min 



Figure 7. Altitude profiles of electron and ion (southern 
streams) density at 20 min in the southern hemisphere. The 
solid line refers to electron density, the dotted line refers to H + 
density, and the dashed line refers to 0 + density. 



-60 -40 -20 0 20 40 60 

latitude (degrees) 

Figure 8. (a) Velocity and (b) density profiles as a function 
of magnetic latitude for the H + southern stream. The solid line 
shows results at 30 min, and the dotted line gives results at 40 
min. 

also moves upward. What is especially interesting, however, 
is that a shock develops near the equator. The steepening be- 
gins at 55 min and is almost complete at 60 min. This shock 
forms as supersonic flow from the southern end is slowed down 
in the northern hemisphere by pressure gradients as well as the 



Figure 9. Density profiles as a function of magnetic latitude 
(between -60 and -40 degrees) at (a) 40 min and (b) 45 min. The 
solid line denotes electrons, the dotted line shows the H + 
southern stream, and the dashed line shows the H + northern 
stream. 
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Figure 10. (a) Velocity and (b) density profiles as a function 
of magnetic latitude for the H+ southern stream. The solid line 
shows results at 50 min, the dotted line gives results at 55 
min, and the dashed line gives results at 60 min. 


electric field. In addition, the flow is accelerated in the south- 
ern hemisphere (between -20 and 0 degrees latitude), primarily 
by the electric field. This shock subsequently moves toward 
the southern end. A similar shock develops in the stream 
which originates in the northern hemisphere. As these shocks 
develop, the electron density between them decreases and con- 
tinues to decrease as they move toward each other. 

The possibility exists that these shocks are artificial, espe- 
cially because in each stream there are large upward flows in 
both hemispheres when they form. Furthermore, when an 
attempt was made to redo the portion of the run from 45 to 60 
min with the cell size equal to about 6.3 km (half the original 
size) and the same Courant number, it was found that very 
narrow, spiky shock structures had developed by 55 min, in 
the approximate place where the large shock formed; when the 
high-altitude time step was reduced, the results were similar to 
those found with the larger cell size. It should be recalled that 
with the larger cell size the large shocks formed at about 60 
min and were not spiky. This enhances the possibility that 
these shocks are numerical artifacts, related to cell size or time 
step problems. Because these shocks are probably artificial, 
further results will not be discussed. 

4. Summary 

We have modeled plasmaspheric refilling following a den- 
sity depletion on an L = 4 field line, using a two-stream ver- 
sion of the plasmaspheric model we developed. The initial 
density profile was chosen to be similar to that of Rasmussen 
and Schunk [1988]. Initially, the upwelling streams pass 
through each other without being impeded. Electrostatic 
shocks are not set up when the streams meet at the equator be- 
cause the stream velocities are highly supersonic with respect 


to the ion-acoustic speed. In addition, Singh [1990] has noted 
that electrostatic shocks are not likely to be set up in hydrody- 
namic models with grid spacings of the order of tens of 
kilometers, since the scale length for the relevant plasma 
processes is of the order of a few Debye lengths (-10 m). It 
should also be noted that perpendicular heating caused by 
wave-particle interactions in the equatorial region is not in- 
cluded. Such heating has been shown [Singh and Torr , 1990; 
Singh and Chan , 1992] to slow down the ion streams enough 
to inhibit the interhemispheric plasma exchange. After the 
stream reaches the conjugate ionosphere, a shock develops 
there which then moves upward. This shock slowly dissipates. 
A reverse shock also forms in each stream in the hemisphere of 
origin and moves up; however, the development of the reverse 
shocks could depend on the initial conditions used. After 1 
hour, large shocks develop in each stream near the equator. 
These shocks then move back toward the equator and down- 
ward. However, these shocks are probably artificial, because 
counterstreaming flows occur in each H + fluid, which the model 
can only handle by creating shocks. 

While the results are similar to those of Rasmussen and 
Schunk [1988], there are differences. Their model is similar to 
our model; however, the ion and electron energy equations are 
included in our model, and our model flux tube extends much 
lower in the ionosphere, to altitudes where 0 + (which is in- 
cluded in our model) is the dominant ion species. The large 
shocks that Rasmussen and Schunk [1988] find are a result of a 
reverse shock which develops in the northern hemisphere, 
moves past the equator, and then reverses direction. Such re- 
verse shocks also develop in our results, but the large shocks 
develop independently of them. After the large shocks form, 
they cross the equator and continue to move downward. 
However, the large shocks that develop in our results are prob- 
ably artificial. 

A similar study was done by Wilson et al [1992], who used 
a semikinetic model that included Coulomb collisions for ion- 
ion interactions. In this model, each simulation particle repre- 
sents many actual particles; the simulation particles are acted 
on by gravitational, magnetic mirror, and electric field forces. 
The electric field is found by using a Boltzmann relation be- 
tween electron number density and electric potential. In the 
model the bases of the flux tube are at about 1900-km altitude. 
At each base, ions are injected into the flux tube, with a veloc- 
ity distribution that corresponds to the upgoing half of a non- 
drifting Maxwellian. Plasmaspheric refilling was studied on a 
variety of L shells, including L = 4. In each case the same se- 
quence of events occurs. Initially, a polar wind type outflow 
from each ionosphere is present. Within an hour, two counter- 
streaming beams exist along the entire flux tube. The phase 
space gap between these two beams is filled in by collisional 
diffusion; when it has filled, the plasma is nearly isotropic and 
Maxwellian everywhere. What is most important, in view of 
the results presented in this paper, is that shock fronts are 
never produced, at the equator or any other point. Field-aligned 
structures in the equatorial region from L = 3 to 5 have, how- 
ever, been inferred from DE 1 observations, especially in the 
electron density and the H + density, temperature, and pitch 
angle distribution [Olsen et al , 1987; Olsen , 1992]. 

Limitations of our model include the restriction to isotropic 
temperatures and the inability to distinguish counterstreaming 
flows within one fluid. In some sense it is inconsistent to 
assume temperature isotropy when using a two-stream model, 
because if the collision frequencies are low enough so that the 
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upwelling streams would not thermalize when they meet, then 
they are low enough for temperature anisotropies to be 
present. The inability to distinguish counterstreaming flows 
within one fluid is a problem if there is a large return flux that 
is about equal in magnitude to the flux flowing into the conju- 
gate hemisphere, since then an artificial shock could be 
formed. Such return fluxes occur after the density of an ion 
stream rises in the conjugate hemisphere and the density gra- 
dient force drives the ions back. The large shocks that form 
near the equator are probably artificial, although this may be 
due to cell size or time step problems rather than the reasons 
given above. However, despite these limitations we have 
shown that initially the two H + streams interpenetrate, and 
that interesting shock structures develop in the ionospheres 
and propagate upward. 
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